#R script GLM regressions

library(brglm2)#v0.9.2; used for Firth regression

GLM_func<-function(metadata_file,PRS,Signature,covariates,output_dir){
  
  #open metafile (PRS+covariates+Signatures)
  metafile<-read_csv(metadata_file)
  
  #run GLM
  result_GLM=data.frame()
  
  for(i in PRS){
    for(j in Signature){
      
      #check the frequency of DNA mutatinal signatures
      eval(parse(text=paste0("len=length(which(metafile$",j,"==1))")))
      
      if(len<=10){
        eval(parse(text=paste0("a<-glm(relevel(as.factor(",j,"),ref='0') ~ scale(",i,")",covariates," , family=binomial ('logit') ,data=metafile)")))
        a0<-update(a, method = "brglmFit")
        a1<-confint.default(a0)
      }
      else{
        eval(parse(text=paste0("a0<-glm(relevel(as.factor(",j,"),ref='0') ~ scale(",i,")",covariates," , family=binomial ('logit') ,data=metafile)")))
        a1<-confint.default(a0)
      }
      
      #Concatenate Results
      rownames(summary(a0)$coefficients)
      temp=data.frame(Variables=rownames(summary(a0)$coefficients),
                      OR=exp(summary(a0)$coefficients[1:length(rownames(summary(a0)$coefficients)), 1]),
                      SE=summary(a0)$coefficients[1:length(rownames(summary(a0)$coefficients)), 2],
                      L95=exp(a1[1:length(rownames(summary(a0)$coefficients)),1]),
                      U95=exp(a1[1:length(rownames(summary(a0)$coefficients)),2]),
                      pval=summary(a0)$coefficients[1:length(rownames(summary(a0)$coefficients)), 4],
                      Signature=j,
                      PRS=i,
                      Model=ifelse(len<=10,'Firth','Ordinary'))
      result_GLM=rbind(result_GLM,temp)
    }
    
  }
  
  #save
  write_csv(result_GLM,paste0(output_dir,'results.csv'))
}

#example of use
#GLM_func(metadata_file ='All_info_PRS.csv',PRS = c('PRS_BMI','PRS_SystolicBP'),Signature = c('SBS1','SBS44'),
#         covariates = '+Age+Sex+PC1+PC1+PC2+PC3+PC4+PC5',output_dir = './')


